OHBM2024 Tutorial: TAC Data Hands-on

Setup

Packages

First we load up the packages we’ll be working with

# remotes::install_github("mathesong/bloodstream")
# remotes::install_github("mathesong/kinfitr")

library(tidyverse)

── Attaching core tidyverse packages ────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.1     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ tidyr     1.3.1
✔ purrr     1.0.2     
── Conflicts ──────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors

library(bloodstream)
library(kinfitr)
library(job)

theme_set(theme_light())

Data

download.file("https://www.dropbox.com/scl/fi/ax1t3lfop3pp3h4cir0m7/ds004869.zip?rlkey=6z7e9wopud9ngu26ekbd3kfks&dl=1", 
              destfile = "ds004869.zip")

unzip(zipfile = "ds004869.zip")

Blood Analysis

job({ 
  bloodstream("ds004869/") 
})

Alternatively, if you want to use a configuration file

bloodstream("ds004869/", configpath = "ds004869/code/config_tutorial.json") 

Loading Data

Loading bloodstream outputs

We want the arterial input function data, which are called “input”.

bloodstream_data <- bloodstream_import_inputfunctions("ds004869/derivatives/bloodstreamtutorial/") %>%
  select(-measurement)

head(bloodstream_data, n = 6)

bloodstream_data %>% 
  slice(1:6) %>% 
  pull(input) %>% 
  map(plot)

[[1]]

Warning: Removed 1 row containing missing values or values outside the scale range (`geom_line()`).

[[2]]

Warning: Removed 1 row containing missing values or values outside the scale range (`geom_line()`).

[[3]]

Warning: Removed 1 row containing missing values or values outside the scale range (`geom_line()`).

[[4]]

Warning: Removed 1 row containing missing values or values outside the scale range (`geom_line()`).

[[5]]

Warning: Removed 1 row containing missing values or values outside the scale range (`geom_line()`).

[[6]]

Warning: Removed 1 row containing missing values or values outside the scale range (`geom_line()`).

Loading petprep outputs (TACs and morphology)

petprep_data <- kinfitr::bids_parse_files("ds004869/derivatives/petprep_extract_tacs/") %>% 
  unnest(filedata) %>% 
  filter(str_detect(path_absolute, "gtmseg"))

tacs <- petprep_data %>% 
  filter(measurement=="tacs") %>% 
  filter(is.na(pvc)) %>% 
  mutate(tacdata = map(path_absolute, ~read_delim(.x, delim="\t", 
                                                  show_col_types = FALSE)))

morphdata <- petprep_data %>%
  filter(measurement=="morph") %>%
  mutate(morphdata = map(path_absolute, ~read_delim(.x, delim="\t",
                                                  show_col_types = FALSE)))

Here, we would usually combine TACs according to which combined regions we were interested in, and use volume-weighted averaging. In this case, we’ll just skip that step and be a little bit naughty and focus on one region: the left putamen.

Creating combined TACs

Now we want to combine regions by volume-weighted averages of the contitutent regions. We’ll create a frontal cortex region, a striatal region and a hippocampus-amygdala region.

First, let’s define the regions

frontal_regions <- morphdata$morphdata[[1]] %>% 
  filter(str_detect(name, "frontal")) %>% 
  pull(name)

frontal_regions

 [1] "ctx-lh-caudalmiddlefrontal"  "ctx-lh-lateralorbitofrontal" "ctx-lh-medialorbitofrontal"  "ctx-lh-rostralmiddlefrontal"
 [5] "ctx-lh-superiorfrontal"      "ctx-lh-frontalpole"          "ctx-rh-caudalmiddlefrontal"  "ctx-rh-lateralorbitofrontal"
 [9] "ctx-rh-medialorbitofrontal"  "ctx-rh-rostralmiddlefrontal" "ctx-rh-superiorfrontal"      "ctx-rh-frontalpole"         

striatal_regions <- morphdata$morphdata[[1]] %>% 
  filter(str_detect(name, "Putamen|Accumbens|Caudate")) %>% 
  pull(name)

striatal_regions

[1] "Left-Caudate"         "Left-Putamen"         "Left-Accumbens-area"  "Right-Caudate"        "Right-Putamen"        "Right-Accumbens-area"

hipamg_regions <- morphdata$morphdata[[1]] %>% 
  filter(str_detect(name, "Hippocampus|Amygdala")) %>% 
  pull(name)

hipamg_regions

[1] "Left-Hippocampus"  "Left-Amygdala"     "Right-Hippocampus" "Right-Amygdala"   

Now we combine the TACs and the region sizes

selected_tacs <- select(tacs, c(ses:rec, tacdata)) %>% 
  inner_join(select(morphdata, c(ses:rec, morphdata))) %>% 
  group_by(sub, ses) %>% 
  mutate(tacdata = map(tacdata, ~pivot_longer(.x, 
                                              cols = `Left-Cerebral-White-Matter`:`ctx-rh-insula`, 
                                              names_to = "name", values_to = "TAC"))) %>% 
  mutate(tacdata = map2(tacdata, morphdata, ~inner_join(.x, .y, by="name")))

Joining with `by = join_by(ses, sub, task, trc, acq, run, rec)`

Then we perform volume-weighted averaging

volume_weighted_average_tac <- function(tacdata, regions, regionname) {
  
  tacdata_combined <- tacdata %>% 
    filter(name %in% regions) %>% 
    group_by(frame_start, frame_end) %>% 
    mutate(volume_tot = sum(`volume-mm3`),
           volume_frac = `volume-mm3` / volume_tot,
           TAC_frac = TAC * volume_frac) %>% 
    summarise(!!regionname := sum(TAC_frac), 
              .groups = "keep") %>% 
    ungroup()
  
  return(tacdata_combined)
  
}

volume_weighted_average_tacs <- function(tacdata, regions_list) {
  
  regionnames <- names(regions_list)
  
  out <- map2(regions_list, regionnames, ~volume_weighted_average_tac(tacdata, .x, .y))
  
  out <- purrr::reduce(out, dplyr::inner_join, by = c("frame_start",
                                                      "frame_end"))
  
  return(out)
}

regions_list <- list(
  Frontal = frontal_regions,
  Striatum = striatal_regions,
  HippAmg = hipamg_regions
)


selected_tacs <- selected_tacs %>% 
  group_by(sub, ses) %>% 
  mutate(selected_tacdata = map(tacdata, ~volume_weighted_average_tacs(.x, 
                                                                       regions_list)))
  

Combining the TAC and blood (input) data

modeldata <- selected_tacs %>% 
  select(ses:rec, tacdata = selected_tacdata) %>% 
  inner_join(bloodstream_data)

Joining with `by = join_by(ses, sub, task, trc, acq, run, rec)`

Fitting TACs

Preparation

Correcting Units

The TAC data are in Bq/mL, and the bloodstream data are in kBq/mL. So we correct this.

modeldata <- modeldata %>% 
  mutate(tacdata = map(tacdata, ~.x %>% 
                         mutate(across(.cols = Frontal:HippAmg, 
                                       ~unit_convert(.x, from_units = "Bq", to_units = "kBq")))))

Adding frame midtimes and durations

modeldata <- modeldata %>% 
  mutate(tacdata = map(tacdata, ~.x %>% 
                         mutate(frame_start = frame_start / 60,
                                frame_end = frame_end / 60) %>% 
                         mutate(frame_dur = frame_end - frame_start,
                                frame_mid = frame_start + 0.5*frame_dur)))

Adding weights

modeldata <- modeldata %>% 
  mutate(tacdata = map(tacdata, ~.x %>% 
                         mutate(meanTAC = rowMeans( .x %>% select(Frontal:HippAmg) )) %>% 
                         mutate(weights = weights_create(t_start = frame_start,
                                                         t_end = frame_end,
                                                         tac = meanTAC))))

Fitting a single TAC

fit_2tc <- twotcm(
       t_tac = modeldata$tacdata[[1]]$frame_mid, 
       tac = modeldata$tacdata[[1]]$Frontal, 
       input = modeldata$input[[1]], 
       weights = modeldata$tacdata[[1]]$weights, 
       vB = 0.05, multstart_iter = 5)

fit_2tc$par

plot(fit_2tc)

Fitting multiple TACs

Here we’ll use a linearised model because they fit more quickly, in this case Logan. But most linearised models require a t* value to operate.

Selecting a t* value

Let’s choose an appropriate t* value

modeldata %>% 
  ungroup() %>% 
  filter(ses=="baseline") %>% 
  slice(1:5) %>% 
  mutate(tstarplot = map2(tacdata, input, 
     ~Logan_tstar(
         t_tac = .x$frame_mid, 
         lowroi =  .x$HippAmg,
         medroi =  .x$Frontal,
         highroi = .x$Striatum, 
         input = .y,
         vB = 0.05)
     )) %>% 
  pull(tstarplot)

Warning: There were 15 warnings in `mutate()`.
The first warning was:
ℹ In argument: `tstarplot = map2(...)`.
Caused by warning:
! Removed 1 row containing missing values or values outside the scale range (`geom_point()`).
ℹ Run ]8;;ide:run:dplyr::last_dplyr_warnings()dplyr::last_dplyr_warnings()]8;; to see the 14 remaining warnings.

[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

Ok, let’s use 13 frames.

Fitting

Let’s focus on the frontal cortex

modeldata <- modeldata %>% 
  group_by(sub, ses) %>% 
  mutate(Logan_fit = map2(tacdata, input, ~Loganplot(t_tac = .x$frame_mid,
                                             tac = .x$Frontal, 
                                             input= .y, 
                                             tstarIncludedFrames = 13, 
                                             weights = .x$weights, 
                                             vB = 0.05, 
                                             dur = .x$frame_dur)))

Plotting

Let’s see a few fits

map(modeldata[1:6,]$Logan_fit, plot)

[[1]]

Warning: Removed 1 row containing missing values or values outside the scale range (`geom_point()`).

[[2]]

Warning: Removed 2 rows containing missing values or values outside the scale range (`geom_point()`).

[[3]]

Warning: Removed 2 rows containing missing values or values outside the scale range (`geom_point()`).

[[4]]

Warning: Removed 1 row containing missing values or values outside the scale range (`geom_point()`).

[[5]]

Warning: Removed 1 row containing missing values or values outside the scale range (`geom_point()`).

[[6]]

Warning: Removed 1 row containing missing values or values outside the scale range (`geom_point()`).

Outcomes

Logan_outcomes <- modeldata %>% 
  select(sub, ses, Logan_fit) %>% 
  mutate(Vt = map_dbl(Logan_fit, c("par", "Vt"))) %>% 
  select(-Logan_fit)

ggplot(Logan_outcomes, aes(x=ses, y=Vt, colour=sub, group=sub)) +
  geom_point(size=3, colour="black") +
  geom_point(size=2.5) +
  geom_line() + 
  scale_y_log10()